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Beschreibung 

Rechnergestiitztes Simulationsverf ahren zur Bestinnnung eines 
elektromagnetischen Feldes eines Korpers 

Die Erfindung betrifft ein rechnergestiitztes Simulations- 
verfahren zur Bestinnnung eines elektromagnetischen Feldes 
eines Korpers, der mehrere Teilgebiete umfaSt und mehrere 
Ladungen und Strome enthalt . 



Aus der DIN-Norm VDE 0 870 ist die "Elektromagnetische 
Vertraglichkeit" (EMV) einer elektrischen Einrichtung bekannt 
als die "Fahigkeit einer elektrischen Einrichtung, in ihrer 
elektromagnetischen Umgebung zuf riedenstellend zu 
15 f unktionieren, ohne diese Umgebung, zu der auch andere 
Einrichtungen gehoren, unzulassig zu beeinf lussen" . 

Ein zu der Klasse der CG-Verf ahren (=Verf ahren der 
konjugierten Gradienten, engl . : conjugate gradient) 
2 0 gehorender Algorithmus mit der Bezeichnung GMRES-Verf ahren 
ist aus [3] bekannt. 

Die schnelle Multipolmethode, nachfolgend als FMM (engl.: 
fast multipole method) bezeichnet ist z.B. aus [4] bekannt. 
Ein wesentlicher Nachteil bei diesem Formalismus ergibt sich 
aus der Eigenschaft, daS die Multipolkoef f izienten nicht 
explizit berechnet werden, wodurch storende Verzerrungen 
auftreten konnen, die wiederum zu Fehlern in den 
Multipolentwicklungen fuhren. 

30 

Die Auf gabe der Erfindung besteht darin, das 
elektromagnetische Feld eines Korpers zu bestimmen durch 
einen Rechner, um bereits vor der Realisierung 
elektromagnetische Vertraglichkeit des Korpers gewahrleisten 
35 zu konnen und somit aus Messungen und Nachbesserungen 

bestehende Opt imie rungs zykl en zu vermeiden, wobei die oben 
dargestellten Nachteile vermieden werden. 
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Die Aufgabe der Erfindung wird gemaS den Merkmalen des 
Patentanspruchs 1 gelost . 

5 In einem rechnergestiitzten Simulationsverf ahren wird ein 
elektromagnetisches Feld eines Korpers, der mehrere 
Teilgebiete umfaSt und mehrere Ladungen und Strome enthalt, 
bestimmt . Dazu werden in den Teilgebieten eine globalen 
Multipolentwicklung und eine lokale Multipolentwicklung 

10 durchgef iihrt , wobei eine globale Multipolentwicklung die 

Wirkung der Ladungen und Strome fur entfernte Punkte in Form 
einer Multipolentwicklung darstellt und analog dazu eine 
lokale Multipolentwicklung die Wirkung der Ladungen und 
Strome an Punkten innerhalb eines der mehreren Teilgebiete in 

15 Form einer Multipolentwicklung darstellt. Globale 

Multipolkoef f izienten werden als Parameter der globalen 
Multipolentwicklung berech.net, lokale Multipolkoef f izienten 
werden als Parameter der lokalen Multipolentwicklung 
berechnet. Schliefilich wird das elektromagnetische Feld des 

2 0 Korpers aus der globalen und der lokalen Multipolentwicklung 
bestimmt . 

Mit der Simulation von fur EMV notwendige vorgeschriebene 
Messungen mittels numerischer Verfahren wird die Optimierung 
25 der Produktkomponenten wesentlich erleichtert . Der sich 

wiederholende ProzeS von Messung und Nachbesserung, wobei die 
simulierte Messung an alien Orten moglich ist, und damit 
verbundene Kosten und Zeit lassen sich mit dem 
erf indungsgemaSen Verfahren wesentlich reduzieren. 

30 

Das hier vorgestellte Verfahren eignet sich insbesondere zur 
Simulation von Streuproblemen, bei denen die geometrischen 
Strukturen der Modelle und damit auch die Dimensionen der 
Teilgebiete kleiner als die Wellenlange sind. Unter 
35 derartigen Voraussetzungen bleiben die Interf erenzef f ekte 
begrenzt, so daS bereits wenige Multipolkoef f izienten zur 
Beschreibung von Potentialen bzw. Feldstarken ausreichen. Mit 
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dem hier beschriebenen Stabilisierungsverf ahren verha.lt sich 
das Verf ahren auch fur tief e Frequenzen stabil . 

Eine erste Weiterbildung der Erfindung besteht darin, den 
5 Korper in mindestens ein weiteres Teilgebiet zu unterteilen. 
Somit besteht der Korper aus Teilgebieten . 

Eine andere Weiterbildung besteht darin, iterativ mit einer 

Startbedingung 1=0 bis zu einer Abbruchbedingung, die ein 
10 ausreichend kleines Fehlermafi vorgibt, folgende Schritte zur 

Bestimmung des elektromagnetischen Feldes des Korpers 
^| durchzuf iihren : 

a) Berechnung der globalen Multipolentwicklungen mit globalen 



Multipolkoef f izienten gemaS 



15 



GI , 



wobei 




20 



c u ein Vektor ist, der die globalen 

Multipolkoef f izienten der mehreren Teilgebiete 
umf aSt , 

I ein Vektor ist, der eine gegebene Stromverteilung 

angibt, 

G eine Matrix ist, die eine Bestimmung der globalen 



Multipolkoef f izienten in dem einen Teilgebiet der 
mehreren Teilgebiete fur die gegebene 
Stromverteilung I angibt; 



b) Berechnung der lokalen Multipolentwicklung mit lokalen 
Multipolkoef f izienten gemafi 



30 




Tc9, 



35 



wobei 
c 1 



ein Vektor ist, der die lokalen 

Multipolkoef f izienten der mehreren Teilgebiete 

umf aSt , 
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T eine Translationsmatrix ist, durch die die 

globalen Multipole in lokale Multipole 
zusammengef iihrt werden; 
c) Bestimmung des elektromagnetischen Felds aus 



5 



ZI 



= Z 1 I + Lc 



1 



wobei 



10 



L 



Z 



eine Impedanzmatrix, 

einen Teil der Impedanzmatrix Z, der Kopplungen 
zwischen den Teilgebieten darstellt, 
eine Matrix zur Auswertung der lokalen 
Multipolkoef f izienten, 



bezeichnen . 



15 



Eine nachste Weiterbildung besteht darin, die Teilgebiete 
gleich grofi zu wahlen. 

In einer Weiterbildung wird die GroSe der Teilgebiete 
20 proportional zu einer Entfernung von einem Beobachtergebiet 



Auch eine Weiterbildung liegt in der Zuordnung jeweils eines 
Teilgebiets zu einem Bereich mit einheitlicher physikalischer 
25 Eigenheit. Dies konnen groSe Flachen mit gleicher 

elektromagnetischer Eigenschaft, beispielsweise gleicher 
Ladungsverteilung o.a., sein. 

Im Rahraen einer anderen Weiterbildung werden 
30 impedanzbehaf tete Bauteile innerhalb des zu untersuchenden 
Korpers direkt als eine Impedanz in der Matrix Z' 
berucksichtigt . 

Ferner kann das elektromagnetische Feld fur vorgebbare 
35 Frequenzen bestimmt werden. 



gewahlt . 
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Eine zusatzliche Weiterbildung besteht darin, die vorgebbaren 
Frequenzen von einer Jcleinsten zu untersuchenden Frequenz bis 
hin zu einer groSten zu untersuchenden Frequenz fortgesetzt 
mit einer vorgebbaren Schrittweite zur Bestiraraung des 
5 elektromagnetischen Feldes auszuwahlen. Dabei kann bei der 
kleinsten zu untersuchenden Frequenz begonnen werden und bis 
zur grofiten zu untersuchenden Frequenz fortgefahren werden 
oder umgekehrt . Zusatzlich kann auch zwischen der kleinsten 
und der groSten zu untersuchenden Frequenz begonnen werden, 
10 wobei bis hin zur kleinsten oder zur groSten Frequenz das 
Verfahren zur Bestimmung des elektromagnetischen Feldes 
fortgesetzt wird. 



Eine zusatzliche Weiterbildung besteht darin, Stabilitat des 
15 Verfahrens bei tiefen Frequenzen dadurch zu gewahrleisten, 
dafi die globalen Multipolentwicklungen elementweise anstatt 
nach Basisf unktionen durchgefuhrt wird. 

Eine Weiterbildung besteht darin, das Verfahren zur 
20 Bestimmung einer elektromagnetischen Vertraglichkeit eines 
Korpers zu verwenden . 



4 



Weiterbildungen der Erfindung ergeben sich auch aus den 
abhangigen Anspriichen. 

Anhand der folgenden Figuren werden Ausfiihrungsbeispiele der 
Erfindung naher dargestellt. 



Es zeigen 

30 Fig.i Ein Blockdiagramm, das Schritte eines Verfahrens zur 
Bestimmung des elektromagnetischen Feldes eines 
Korpers zeigt, 

Fig. 2 einen ideal leitenden Korper in einem auSeren Feld, 
Fig. 3 geometrische GroSen zu einem Linien- 
35 Gewichtungsverf ahren, 

Fig. 4 einen EinfluS eines n-ten Elementpaars auf ein m-tes 
Elementpaar , 
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Fig. 5a eine verteilte Impedanz bei impedanzbehaf teten 
Elementen, 

Fig. 5b eine konzentrierte Impedanz bei impedanzbehaf teten 
Elementen, 

5 Fig. 6 eine Einteilung einer Autokarosserie in Teilgebiete 
fur die Verf einerungsstuf en 2,3 und 4, 
Fig. 7 eine Einteilung in Teilgebiete bei variierenden 
Element grofien, 

Fig. 8 eine Skizze, die Kugelkoordinaten r, S, a darstellt, 
10 Fig. 9a eine Skizze, die geometrische Grofien zur globalen 
Multipolentwicklung darstellt, 
Fig. 9b eine Skizze, die geometrische Grofien zur lokalen 

Multipolentwicklung darstellt, 
Fig. 10a eine Skizze, die eine Leiterschleif e und eine 
15 Einteilung der Basisfunktionen in vier Teilgebiete 

darstellt , 

Fig. 10b eine Skizze, die an den Randern der Teilgebiete 



Fig. 11 eine Skizze, die direkt und indirekt gekoppelte 

Teilgebiete darstellt, 
Fig. 12 eine Skizze, die geometrische Grofien zur Translation 

darstellt, 

25 Fig. 13 eine Skizze, die direkt und indirekt gekoppelte 
Teilgebiete fur eine zweistufige FMM darstellt. 

Zunachst wird die Erfindung anhand der Figuren 2 bis 13 naher 
erlautert . 

30 

DIE RETARD I ERTEN POTENTIALE 

Ausgegangen wird von den Maxwellschen Gleichungen in 
dif f erentieller Form, wie sie dem Fachmann hinreichend 
bekannt sind: 



20 



entstehende Ladungen, die sich nicht mehr ausloschen 
und zu Abbruchf ehlern in der Entwicklung von cp 
fiihren, zeigt, 



35 
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V • E = — (la) 

e 0 

V x E = -i©Ho H 

5 V • H = 0 (lc) 

V x H = J + ia>s 0 E (Id) , 
wobei 

10 V ein Nabla-Operator , 

H eine magnetische Feldstarke, 

E eine elektrische Feldstarke, 

p eine Raumladungsdichte, 

Sq eine Permittivitat des Vakuums, 

15 i eine imaginare Einheit (V-l), 

co eine Kreisf requenz , 

MO eine Permeabilitat des Vakuums, 

J eine Volumenstromdichte 

bezeichnen . 

20 

Eine Kontinuitatsgleichung folgt unmittelbar aus den 
Maxwellschen Gleichungen durch Kombination von Gleichung (la) 
mit der Divergenz von Gleichung (Id) : 

25 V • J = -icop (2) 

Werden die Potentialansatze 

E = -icoA - Vcp (3a) 



30 



H = — V x A (3b) 
HO 

verwendet, ergeben sich unter Berucksichtigung der Lorentz- 
Ei chung 



35 
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V • A + icoeQUo^ = 0 ^ 
die beiden Dif f erentialgleichungen : 

(v 2 + k 2 ^A = -n 0 J (5a) 



und 



(v 2 + k 2 ) q> = - — (5b) 
V I e 0 

10 

mit der Wellenzahl 

I CO 271 , , 

k = co^/eo^O = — = — (6) , 

c 0 A, 

15 wobei 

A ein magnetisch.es Vektorpotential , 

<p ein elektrisches Skalarpotential , 

Cq eine Vakuumlichtgeschwindigkeit , 

X eine Wellenlange 

20 bezeichnen. 

Sowohl das skalare elektrische Potential q> als auch die drei 
kartesischen Komponenten des magnetischen Vektorpotentials A 
werden durch dieselbe Art von Dif f erentialgleichung 
25 beschrieben, die i.a. als Helmholtz-Gleichung bezeichnet 
wird : 

(v 2 + k 2 ) M/(r) = f(r) (7) , 

30 wobei 

y/ eine Basisf unktion, 

r einen Ortsvektor, 

f () eine Funktion von einem Argument 

bezeichnen. 
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25 



1st die Losung der Helmholtz -Gleichung fur eine diracformige 
Anregung auf der rechten Seite bekannt, so kann die 
Gesamt losung fiir eine gegebene Funktion f (r) durch 
Superposition bestimmt werden. Aus Symmetriegrunden ergibt 
sich fiir eine diracformige Anregung im Ursprung eine 
kugelsymnietrische Losungs funktion if/, so daE die Gleichung 
(7) in Kugelkoordinaten auf die von r abhangige 
Dif f erentialgleichung 

1 d 1 d 77 

— — \|/ + kV* = 0 (8) 
r ^ or ox 



reduziert wird. Wird in Gleichung (8) die Substitution 
15 vj> = rv|/ 

durchgefiihrt , so folgt eine gewohnliche Dif f erentialgleichung 
mit konstanten Koef f izienten : 

2 0 vj>" + k 2 vp = 0 (9) 

mit der allgemeinen Losung: 

vj> = Cje -11 ^ + C 2 e lkr (10) 



beziehungsweise 



e ~ikr e ^" kr 
M/ = C x + C 2 (11) 



30 Der erste Term aus Gleichung (11) beschreibt eine 

auslaufende, der zweite Term eine einlaufende Welle. Da 
letztere aus physikalischen Griinden ausscheidet, reduziert 
sich das weitere Vorgehen auf die Bestimmung der Konstanten 
C-l . Diese ergibt sich fiir das elektrische Skalarpotential (p 
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durch Vergleich der allgemeinen Losung mit dem 
elektrostatischen Potential (k— >0) eines sich am Ursprung 
befindlichen Ladungselements pdV zu: 

5 Cf 1 = — ^- pdV (12) . 

x 4tu£ 0 

Analog ergibt sich die Konstante fur das magnetische 
Vektorpotential eines sich im Ursprung befindlichen 
Stromelements JdV zu : 

0 

c mag = iiQ. Jdv (13) 

^ 471 



Durch Integration iiber samtliche Quellgebiete G ergeben sich 
15 die Losungen der Gleichungen (5a) und (5b) schlieSlich als: 

— ikR 

*W " ^JJJ^O^-^-IV ,14a) 
G 

<p(r) = fffp( r ') dV (14b), 

U G 



20 



mit : R = fir - r' 



wobei 

R einen Abstand zwischen Beobachterpunkt r und 

25 Quellpunkt r' , 

G ein Volumengebiet , 

A, <p retardierte Potentiale 

bezeichnen. 
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BERECHNUNG PER FELDSTARKEN 

Durch Einsetzen der Gleichung (14a) in (3a) und der Gleichung 
(14b) in (3b) lassen sich die elektrische Feldstarke und die 
magnetische Feldstarke bestimmen: 

E( r ) = _ fff J(r , } e dv , _ _±_ fff (r , )v a dv , 

471 JJJ R 4718 0 JJJ R 

G G 

(15a) 

H(r) = — fff V x j(r') dV (15b) 

47t JJJ R ' 

G 

Unter Verwendung von Kugelkoordinaten f olgt somit : 



„ e ikR -ikRf 1 ik'. 

V — = -e + -jlR (16) 



Die Integralf ormeln zur Bestimmung der elektrischen und 
magnetischen Feldstarke bei gegebenen Quellen ergeben sich 
zu : 

^ - - ^ f ™ v- *" + ±- 0 f ^ ^ + S r 

(17a) 

H(r) = ± fff e- ikR f-^- + i±) J(r-) x R dV (17b) 
471 „ ^R R ^ 



Die Gleichung (17b) kann als Erweiterung des Biot-Savartschen 
Gesetzes fur den zeitharmonischen Fall aufgefafit werden. 
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IDEALE LETTER IM AUSEREN MAGNET I S CHEN FELD 

Wird ein ideal leitender Korper in ein auEeres einfallendes 

Feld E 1 , H 1 gebracht, so bildet sich auf der Oberf lache S des 
ideal leitenden Korpers eine Strom- und Ladungsverteilung 
K, a derart, dafi das Innere des Korpers feldfrei wird (siehe 
Fig. 2 ) . Das von den induzierten Quellen K, a erzeugte Feld 

wird als Streufeld E s , H s bezeichnet . Es la£t sich mit Hilfe 
der Gleichungen (17a) und (17b) berechnen und ist dem 

einfallenden Feld E 1 , H 1 iiberlagert . Da das Gesamtfeld 

E = E 1 + E S (18) 

H = H 1 + H S (19) 

15 im Inneren des idealen Leiters verschwindet , ergeben sich die 

Randbedingungen auf der Innenseite der Leiteroberf lache S 
zu : 



10 



20 



25 



E tan = 0 fur r e S (2 0a) 

H tan = 0 fur r e S~ (20b) . 

Beim Durchlaufen der Randf lache verhalt sich H tan unstetig, 

so daS fur die Aufienseite der Leiteroberf lache S + gilt: 

E tan = 0 fur r e S + (21a) 

H tan = K x n fur r e S + (21b) . 

30 Dabei ist die Flachenstromdichte K mit der 

Flachenladungsdichte a liber die zweidimensionale 
Kontinuitatsgleichung 



V s • K = 



-icoa 



(22) 
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verkniipf t . Der Ausdruck Vg • K bezeichnet die 
zweidimensionale Divergenz der Flachenstromdichte K in der 
Randf lache S . 

5 

Befinden sich mehrere ideal leitende Korper im einfallenden 
Feld, so setzt sich die Randflache S aus verschiedenen 
Teilflachen Sj^ zusammen, also 

10 S = (J S k (23) . 



Das Streufeld E s , H s umf aSt dann samtliche Teilfelder, die 



von den einzelnen Korpern ausgehen, die Randbedingungen gemaS 
den Gleichungen (20a) , (20b) , (21a) und (21b) gelten fur alle 
15 Leiteroberf lachen . 

Die Aufgabe eines Streuproblems besteht nun darin, fur 

gegebene einf allende Felder E 1 , H 1 unter Einhaltung der 
Randbedingungen (20a, b), (2la,b) die Strom- und 
20 Ladungsverteilung K, a auf den Oberflachen des Leiters zu 
f inden . 

Nachfolgend wird der Begriff "Stromverteilung" als 
Bezeichnung fur eine gekoppelte Strom- und Ladungsverteilung 
25 benutzt, ohne explizit auf die Ladungsverteilung einzugehen. 
Diese laSt sich jederzeit mit Gleichung (22) herleiten und 
kann somit dem Begriff "Stromverteilung" untergeordnet 
werden . 

30 Fur das Feld auBerhalb der ideal leitenden Korper kann ein 
klassischer Eindeutigkeitssatz herangezogen werden: 

Die Feldverteilung in einem geschlossenen Volumengebiet G 
ist eindeutig durch die darin enthaltenen Quellen und 
35 durch den Verlauf von E tan oder H tan auf dem Rand des 

Volumengebiet s G spezif iziert . 
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Fur die Eindeutigkeit behandelter Streupr obi erne reicht die 
Einhaltung einer der Randbedingungen der Gleichungen (20a, b) 
bzw. (21a, b) aus . 

Eine Auswertung von E tan durch Kotnbination der Gleichung 
(18) mit der Gleichung (20a) ergibt 

-n x E s = n x E 1 fur r e S (24) 

Da sich E tan beim Durchlaufen der Randflachen stetig 
verhalt, spiel t es hierbei keine Rolle, ob die Beziehung auf 

der AuEenseite der Leiteroberf lache S + oder auf der 
Innenseite der Leiteroberf lache S - ausgewertet wird. Das 
15 Streufeld E s laSt sich mit Gleichung (17a) bestimmen. Wird 
weiterhin die Flachenladungsdichte a, unter Verwendung der 
Kontinuitatsgleichung (22) , durch die Flachenstromdichte K 
ausgedriickt, so folgt die elektrische Feldintegralgleichung 
in integrodif f erentieller Form: 



10 



20 



e -ikR 

n x 



471 
S 



da' = n x E^r) 



icou 0 K(r') + V s • K(r') ( _i_ + ik" 
R icos 0 l R 3 R 2 y 

fur reS (25) . 



Eine Auswertung der tangentialen magnet ischen Feldstarke H tan 
25 durch Kotnbination der Gleichung (19) mit Gleichung (20b) bzw. 
Gleichung (21b) fuhrt zu den Beziehungen 

-n x h s = n x H 1 fur r e S~ (26a) 

30 K - n x H s = n x H 1 fiir r e S + (26b) 



Anders als bei Gleichung (24) muS hier berucksichtigt werden, 
auf welcher Seite der Flache S sich der Beobachterpunkt r 
befindet. Urn zu einer davon unabhangigen Integralgleichung zu 
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gelangen, mu£ das Sprungverhalten von H tan beim Durchlaufen 
der Flache S naher betrachtet werden. 

Es folgt schlieSlich die magnetische Feldintegralgleichung 
fur ideal leitende Korper mit glatten Oberflachen: 

^K(r) - n x Jj£±^U- + i|) K(r') x R da- = n x H^r) 

2 J J 471 R <*J 

fur r e S (27) 



Da der Einsatz der magnet ischen Feldintegralgleichung (27) 
sich i.a. auf Korper mit glatten Oberflachen und grofien 
Krummungsradien beschrankt, wird in der vorliegenden 
Erfindung weiterhin von der elektrischen 

Feldintegralgleichung ausgegangen. Dabei wird nicht direkt 
auf Gleichung (25) zuriickgegrif f en, sondern eine leicht 
modifizierte Formulierung der Gleichung (25) verwendet: 



n x 



471 



— K(r,) iss^ V W — v * ■ K <"> * 



fiir r e S 



= n x E-'-(r) 

(28) 



Die Losung der vier Maxwellschen Gleichungen unter Beachtung 
der Randbedingungen (20a, b), (21a, b) reduziert sich somit auf 
die Analyse der Integralgleichung (28). Die dabei gesuchte. 
Losungsf unktion ist die Stromverteilung K(r) auf den 
Oberflachen der beteiligten Korper. 



VEKTORBAS I S FUNKT I ONEN 

Fiir die Stromverteilung K, a ergibt sich unter 
Berucksichtigung der Kontinuitatsgleichung (22) : 
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N 



K(r) = a nM>n( r ) 



(29a) 



n = l 



1 



N 



a(r) 



ico 



Z a n V S * M>n(*) 



(29b) , 



n = l 



5 



wobei 



10 



N 



n 



a 




n 



eine Vektorbasisf unktion, 

einen zweidimensionalen Nabla-Operator in der 
Flache S, 

ein Moment (skalarer Koeffizient) , 
eine Anzahl der Basisf unktionen, 
eine Zahlvariable 



bezeichnen . 

Die Verwendung von Vektorbasisf unktionen auf Dreieckelementen 
15 ist aus [1] bekannt . In [2] sind Vektorbasisf unktionen fur 
parametrische Elemente dargestellt. Darin ist die Idee 
enthalten, mit Hilfe einer linearen Flachenstromverteilung 
eine stiickweise konstante Ladungsverteilung zu erreichen. Da 
sich Vektorbasisf unktionen erster Ordnung iiber zwei 
20 benachbarte Flachenelemente erstrecken, wirken pro 

Zylinderelement bis zu zwei, pro Dreieckelement bis zu drei 
und pro Rechteckelement bis zu vier Basisf unktionen zusammen. 
Eine Regel fur die Zuordnung von Basisf unktionen zu 
Flachenelementen lautet demnach: 



30 Vektorbasisf unktionen erster Ordnung verfugen i.d.R. uber 
folgende Eigenschaf ten : 

a) Die sich ergebenden Flachenstromdichte K ist stetig. 

b) Die sich ergebende Flachenladungsdichte a ist 
stiickweise konstant . 



25 



Jeder inneren Kante eines Diskretisierungsnetzes 
entspricht eine Basisf unktion und somit ein Freiheitsgrad 
des resultierenden Gleichungssystems . 



5 
15 

4 

25 
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c) Auf dem Rand des Def initionsgebiets einer 

Basisf unktion gilt Kj_=0 

d) Jede Basisf unktion erfiillt fur sich die 
Kontinuitatsgleichung (22) . 

e) Es treten keine singularen Ladungen auf, wie sie bei 
Basisf unktionen nullter Ordnung unvermeidbar sind. 

f) Jede Basisf unktion kann in groSer Entfernung als ein 
Hertzscher Dipol angesehen werden. 



LINIEN-GEWICHTUNGSVERFAHREN 

Bei dem Linien-Gewichtungsverf ahren wird die Randbedingung 
E tan =0 in Form eines Linienintegrals 

J E • dr = 0 (30) 

vorgegeben. Die Kurve C m bezieht sich dabei auf das m-te 
Elementpaar, bestehend aus den Elementen S ml und S m2 , und 

verlauft geradlinig vom Schwerpunkt des ersten Elements 
zum Mittelpunkt der gemeinsamen Kante und anschlieSend 

^m2 



geradlinig weiter zum Schwerpunkt des zweiten Elements r5 



(siehe Fig. 3 ) . 

Die Idee beim Linien-Gewichtungsverf ahren besteht darin, die 
Gradientenbildung des Skalarpotentials cp durch ein 
Linienintegral in eine diskrete Potentialdif f erenz zu 
iiberfuhren. Ausgegangen wird von der elektrischen 
Feldintegralgleichung (28) 

(icoA + Vcp) tan = E^ an fur r e S (31) 



mit den Potentialen 
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a = Jia 



-ikR 



47t JJ R 
S 



K(r') da' 
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(32a) 



» - --Vlf 

*i mc «. Air J J 



-ikR 



10Sn47l JJ R 
U S 



V s • K(r-) da' 



(32b) 



5 Werden die beiden Seiten aus Gleichung (31) vektoriell iiber 
die M Kurven gewichtet, so ergeben sich die Gleichungen 



| (iooA + Vq>) • dr = J E 1 • dr mit m=l, 2, . . . , 



M 



(33) 



-m 



-m 



10 Die dabei auftretende Integration bezuglich Vq> kann nun 

durch die Potentialdif f erenz der beiden Endpunkte ersetzt 

i 

werden. Werden E und A weiterhin in erster Naherung als 

konstant innerhalb eines Flachenelements angesehen, so ergibt 
sich fur Dreieckelemente : 



15 



^ - 4-) Y + 4&) - - 



mit m=l , 2 , . . . , M 



l m2 



(34) 



2 0 Die beiden lokalen Ortsvektoren dJL und sind in Fig. 3 



dargestellt . 



25 



Ein Vorteil des Linien-Gewichtungsverf ahrens besteht in der 
Umgehung der starken Singularitaten aus Gleichung (25) , 
wodurch sich die Berechnurig erheblich vereinf acht . 
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AUFSTELLEN PES LINEAREN GLE I CHUNGS SYSTEMS 

Fur die Potentiale aus Gleichung (32a) und Gleichung (32b) 
folgen unter Verwendung der Gleichung (29a) und der Gleichung 
(29b) die Beziehungen: 

N -ikR 

A(r) = ijfl. 2 «n JJ VnM da' (3 5a) 

n = 1 ^ 

1 n e - ik R 

q>(r) = — Z a njf— " V s • M> n (r') da' (3 5b) 

n = l Sn 

mit S n : = Sm ^ Sn2 . 

Somit und unter Beriicksichtigung der Beziehung l n = a n lafit 
sich die Gleichung (34) als Matrixgleichung der Form 

Z • I = U (36) " 

darstellen. Dabei bezeichnet die Matrix Z = [Zmn] eine 
Impedanzmatrix, da sie die Stromamplituden I = [l n ] mit den 

Werten U = [u m ] verkniipft. Analog wird die inverse Matrix Z -1 
als Admittanzmatrix bezeichnet. 

Die gesuchten Momente I n lassen sich durch direktes oder 
iteratives Losen des Gleichungs systems (36) bestimmen. 

Nach Aufspaltung der Flache S n in die beiden Elementf lachen 
S nl und S n2 folgt aus Gleichung (34) fur das Matrixelement 



J mn 



Zmn = 



(A 11 + A 21 )4 1 -(Al2+A 22 )^2 



2 v ' 2 

+ <Pl2 + <P22 - <Pll - <P21 (37) 
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wobei 



A 



pq 




vj/ n (r') da' 



(38a) 



5 



<Ppq 



i(OSo47C 



1 




V S • M> n (r') da' 



(38b) 



die einzelnen Potentialbeitrage bezeichnen. Das Matrixelement 
Z,^, auch Koppelimpedanz genannt, beschreibt den EinfluS des 
n-ten auf das m-te Elementpaar . In Fig. 4 sind Koppelpfade 
10 dargestellt zwischen einem Quellelement n und einem 



S n2 erzeugte Vektorpotential im Mittelpunkt von S ml . 

Ein programmtechnischer Ansatz zur Berechnung der 
15 Impedanzmatrix [Z tnn ] besteht in einer separaten, 

vollstandigen Auswertung der einzelnen Terme gemaS 

Gleichung (37) . Da im Fall von Dreieckelementen jedoch bis zu 
drei Basisf unkt ionen pro Element beteiligt sind, liegt es 
nahe, die Berechnung der Impedanzmatrix [Z,^] stufenweise 
20 durchzuf iihren . Dazu werden Teilergebnisse gemaS Gleichung 

(38a) und Gleichung (38b) jeweils fur zwei Dreieckselemente 
berechnet und, mit geeigneten Faktoren versehen, auf die 
betroffenen Matrixelemente addiert. 

25 Die rechte Seite des Gleichungssystems (36) wird gemaS 
Gleichung (34) wie folgt berechnet: 



Beobachterelement m. Beispielsweise bezeichnet A 2 i das von 




(39) . 



1 



5 



15 



20 



25 
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Fur eine einfallende ebene Welle mit k als Ausbreitungsvektor 

i 

und Eq als Polarisationsrichtung lafit sich E fiber die 
Beziehung 

E x (r) = E 0 • e ~ lk ' r (40) 
bestinunen . 



IMPEDANZBEHAFTETE ELEMENTE 

Werden diinne Schichten mit Hilfe einfacher Flachen 
modelliert, so lassen sich auch nichtideale Leiter mit 
ohmschen, induktiven Oder kapazitiven Flachenbelagen relativ 
einfach in das vorhandene Gleichungssystem mit einbeziehen. 
Fur derartige Leiter andert sich die einzuhaltende 
Randbedingung von E tan = 0 auf 



E tan = z S • K (41) 

Aufgrund des Skineffekts ist die Flachenimpedanz z s i.a. 
f requenzabhangig. Fur Schichten, deren Dicke 1^ bei der 
gegebenen Frequenz jedoch klein gegeniiber der Skintiefe ist, 
kann die Stromverteilung uber dem Querschnitt als konstant 
angenommen werden. Beispielsweise laSt sich dadurch die 
Flachenimpedanz z s einer widerstandsbehaf teten metallischen 

Schicht bei tiefen Frequenzen naherungsweise iiber die 
Beziehung 




(42) 



berechnen, wobei r sp den spezifischen Widerstand des 
verwendeten Metalls bezeichnet. 
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Durch die neue Randbedingung (siehe Gleichung (41) ) andert 
sich auch das Gewichtungsintegral aus Gleichung (30) zu 

J E • dr = J z s K • dr (43) , 

c m 

wobei die in Gleichung (43) auf der rechten Seite 
entstehenden GroSen in Form von Korrekturtermen im 
Gleichungssystem zu beriicksichtigen sind. 

10 In Fig. 5a sind verteilte Impedanzen gezeigt, wobei samtliche 
das Flachenelement , hier dargestellt mit der Flachenimpedanz 
z s , betreffende Gewichtungsintegrale beeinfluSt werden. Fur 
die Korrekturterme (also die rechte Seite der Gleichung (43) ) 
gilt 



15 



J z s K • dr = X I n J z s M>n ■ dr (44) , 



so daS pro Dreieckelement in der Regel bis zu neun; pro 
Rechteckelement in der Regel bis zu 16 Matrixeintrage in 
2 0 Gleichung (36) zu korrigieren sind. 



Einfacher ist es, konzentrierte Impedanzen in Form von 
Bauelementen zu beriicksichtigen, wie in Fig . 5b dargestellt 
ist. Im Gegensatz zu verteilten Impedanzen (siehe Fig. 5a) 
25 lassen sich konzentrierte Impedanzen durch eine einzige 

impedanzbehaf tete Kante darstellen. Es wird nur ein einziges 
Gewichtungsintegral und damit auch nur eine einzige Zeile im 
Gleichungssystem (36) beeinf luSt . Fur die mit der Impedanz 

behaftete m-te Kante ergibt sich 

30 

N 

Z z mn I n = U m - (45) , 

n = l 
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wobei die Grofie Z^ in Gleichung (37) und die Grofie U m in 
Gleichung (39) definiert sind. Nachdem der zusatzliche Term 
auf die linke Seite gebracht worden ist, kann dieser in Form 
eines korrigierten Diagonalelements 

^mm = z mm + Zrh (46) 
beriicksichtigt werden. 



KONZ ENTR I ERTE SPANNUNGSOUELLEN 

Analog zum Fall konzentrierter Impedanzen wird auch bei 
konzentrierten Spannungsquellen eine einzige Zeile des 
Gleichungssystems (36) beeinflufit: 

N 

ZZmn^ = U m - U m (47) . 

n = l 

Die zusatzliche Spannungsquelle laSt sich unverandert in 

die recht Seite des Gleichungssytems (36) aufnehmen, wodurch 
sich folgende Beziehung ergibt : 

U m = U m - U m (48) . 

BERUCKSICHTIGUNG VON SYMMETRIEN 

Liegt eine symmetrische Anordnung vor, so laSt sich der 
Aufwand zur Losung des Gleichungssystems erheblich 
reduzieren. Im Fall einer spiegelsymmetrischen Anordnung und 
einer symmetrischen Anregung U m laSt sich die Anzahl der 
Freiheitsgrade um die Halfte reduzieren. Handelt es sich 
hingegen um eine unsymmetrische Anregungen, so wird der 
Anregungsvektor in einen symmetrischen und einen 
unsymmetrischen Teil aufgespalten. 
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Mit dem vorgestellten Verfahren lassen sich auch ideal 
leitende Halbebenen berucksichtigen. 

Im Bereich der EMV- Simulation lassen sich mit Hilfe 
derartiger Uberlegungen beispielsweise die metal lischen Boden 
innerhalb der MeShallen berucksichtigen, ohne da£ eine 
Erhohung der Anzahl der Freiheitsgrade erf olgt . 



10 NUMER I S CHE AUSWERTUNG DER INTEGRALE 

Zur Berechnung der Impedanzmatrix ist es erf orderlich, die 
Integrale aus den Gleichungen (38a) und (3 8b) geeignet 
auszuwerten. Sind Quellpunkt und Beobachterpunkt weit genug 
voneinander entfernt, so geniigt es, die Integrale vollstandig 
15 mit Hilfe numerischer Integrationsf ormeln auszuwerten. In 
Fallen, in denen Quellpunkt und Beobachterpunkt nahe 
beieinander liegen, ergeben sich jedoch Probleme aus der 
Singularitat des Integranden, so daS hierfiir spezielle 
analytische Verfahren zum Einsatz kommen. 



20 



25 



Die in den Gleichungen (38a) und (38b) auftretenden Integrale 
vom Typ 



e -ikR 

ff d da' (49a) 

JJ R 
S 

e -ikR 

ff da' (49b) 

JJ R 
S 



besitzen fur R=0 bzw. r=r' eine schwache Singularitat, die 
abgespalten und analytisch integriert werden kann. Die 
Element f lache, uber die sich die Integration erstreckt, ist 
30 dabei ganz allgemein als S bezeichnet . 
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Fur das Integral aus Gleichung (49a) ergibt sich nach einigen 
Umf ormungen : 

-ikR -ikR _ - 

s s 

, v , 



/ „\ -xkR - 



S 

IG 2 



JJ R 

S 



IG 3 

+ (f - r f ) • JJ ^da' (50a) , 



R 
S 



IG 4 

wobei 
f 

r die freie Ecke des jeweiligen Elements (vorzugsweise 

Dreiecks) , 

¥ einen Vektor, der in die Dreiecksebene projiziert 

wurde (angedeutet durch die Tilde) 
bezeichnen . 



Fur das Integral aus der Gleichung (4 9b) ergibt sich die 
analoge Umf ormung : 

e ~ikR e ~ ikR _ 1 _ 

s s ' s 

IG 2 IG 4 
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Die singularen Integrale IG3 und IG4 lassen sich analytisch 
auswerten. Zur Auswertung der Integrale IG-^ und IG2 werden 
nuraerische Integrationsf ormeln verwendet, da die Integranden 
in der gesamten Elementf lache stetig und beschrankt sind. 

Als Grenzwert fur R— >-0 ergibt sich fur den Integranden aus 
IG 2 nach der Regel von de 1 1 Hospital : 

e -ikR _ x - ik ■ e~ ikR 

lim = lxm = -ik (51) . 

R->0 R R->0 1 

Entsprechend f olgt fur den Integranden aus IG^ : 



lim 
R-»0 



e -ikR _ 1 
(*' - r) 



= 0 (52) 



15 VERFAHREN DER KONJUG I ERTEN GRAD I ENTEN 

Neben der Moglichkeit, die sich ergebenden Gleichungssysteme 
(36) mit direkten Verfahren zu losen, bietet sich aus 
Ef f izienzgriinden der Einsatz von iterativen 
Losungsalgorithmen an. Das zu losende lineare 
2 0 Gleichungssystem aus Gleichung (36) wird allgemein als 

Ax = b mit A e c NxN (53) 

dargestellt, wobei C die Menge der komplexen Zahlen 

25 bezeichnet. Charakteristisch fur iterative Losungsverf ahren 

•. (m) 
1st die generierte Folge von Naherungslosungen x , die als 

Iterierte bezeichnet werden. Konvergiert diese Folge mit 
wachsendem m gegen die exakte Losung x, so kann der 
IterationsprozeS nach Erreichen der gewunschten Genauigkeit 
30 abgebrochen werden. 
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Wahrend direkte Losungsverf ahren einen genau determinierten 
Rechenaufwand besitzen, hangt der Aufwand bei iterativen 
Losungsverf ahren unmittelbar von der Zahl der benotigten 
Iterationsschritte ab und ist im Voraus daher zumeist 
5 unbekannt . 

Zur Losung der Gleichungssysteme benotigen direkte 

Losungsverf ahren einen unmittelbaren Zugriff auf die 

einzelnen Matrixelemente . Bei iterativen Losungsverf ahren 

10 hingegen tritt die Matrix A vorzugsweise nur in Form von 

Matrix- Vektor-Produkten auf und braucht daher nicht explizit 

vorzuliegen. Zum Berechnen der Matrix-Vektor-Produkte geniigt 

der Zugriff auf entsprechende Funktionen. Je nach 

Iterationsverf ahren werden pro Iterationsschritt ein Oder 

15 mehrere Matrix-Vektor-Produkte berechnet , so dafi im Fall von 

2 

voll besetzten Matrizen hierfiir in der Regel 0(N ) 

FlieSkommaoperationen benotigt werden. Wichtig bei der 
Verwendung iterativer Verf ahren ist es somit, die Zahl der 
benotigten Iterationsschritte so klein wie moglich zu halten. 



20 



30 



Zu der Klasse der CG- Verf ahren (=Verf ahren der konjugierten , 
Gradient en, engl . : conjugate gradient) gehort auch das GMRES- 
Verf ahren, das aus [3] bekannt ist. Der dem Verf ahren 
zugrundeliegende Algorithmus basiert auf einer vollstandigen 
QR-Zerlegung des sich ergebenden Losungsraums und ist daher 
auSerst robust. Im der vorliegenden Beschreibung wird das 
GMRES- Verf ahren als Vertreter der CG- Verf ahren verwendet, 
wobei die Erfindung nicht auf die Verwendung dieses 
Verfahrens beschrankt ist. 

Um die Losungsqualitat wahrend des Iterationsprozesses 

beurteilen zu konnen, wird, ausgehend von Gleichung (53), ein 
(m) 

Residuum r 



35 



r( m ) = b - AxH 



(54) 



GR 97 P 1722 



28 

(m) 

als ein FehlermaS fur die Iterxerte x emgefuhrt. Der 

(m) 

tatsachliche Fehler e 



>) = x - *H 



(55) 



10 



ist mangels einer exakten Losung x in der Regel nicht 

(m) 

bekannt . Das Residuum r miSt den Fehler auf der rechten 

Seite, d.h. im Fall der Moment enmethode die Genauigkeit der 
Einhaltung der Randbedingungen nach dem m-ten 
Iterationsschritt . Nach Einsetzen der Gleichung (55) in 
Gleichung (54) ergibt sich unter Berucksichtigung von 
Gleichung (53) der Zusammenhang : 



15 



20 



(56) . 



um einen 



Als Abbruchkriterium fiir den iterationsprozeS wird 

vorzugsweise gefordert, dafi eine Residuumnorm 

-4 

vorgebbaren Faktor e, z.B. 10 , verringert worden ist. 
Falls der Iterationsprozefi mit dem Nullvektor als Anfangswert 
gestartet wird, kann somit die Iteration abgebrochen werden, 
sobald die Ungleichung 



.(■») 



r(°) 



HI 



< 8 



(57) 



25 erfullt ist. Prinzipiell sollte das Abbruchkriterium so 
gewahlt werden, da£ der Fehler durch das Abbrechen der 
Iteration knapp unterhalb des Diskretisierungsf ehlers bleibt 



30 



Die Gesamtzahl der bis zum Erreichen des Abbruchkriteriums 
benotigten Iterationsschritte werden mit M bezeichnet, wobei 
gilt m=l, 2, . . . ,M. 
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VORKONDITIONIERUNG 

Bei den meisten Problemen kann das Konvergenzverhalten des 
Iterationsprozesses durch eine sogenannte Vorkonditionierung 
beschleunigt werden. 

Es kann unterschieden werden zwischen einer 

Vorkonditionierung von links und einer Vorkonditionierung von 
rechts. Aus Gleichung (53) folgt fur das linkstransf ormierte 
Glei chungs system 

( M -1 Ax = M -1 b <=> Ax = b mit A = M _1 A 

b = M -1 b (58) . 

-1 

Die Matrix M wird als Vorkonditionierer bezeichnet . Im Fall 

15 M=A sind beispielsweise samtliche Eigenwerte von A eins, so 
daS sich die exakte Losung bereits nach einer einzigen 

-1 

Iteration einstellt. Der Aufwand fur die Berechnung von M 
ist bei einer derartigen Vorkonditionierung allerdings genau 
so groS wie bei einer direkten Losung des Gl ei chungs syt ems . 
20 Prinzipiell gilt also, dafi je ahnlicher sich die Matrizen A 
und M sind, desto besser ist auch der Effekt der 
Vorkonditionierung . 



10 

* 



# 



Bei Anwendung der Gleichung (58) muS die GroEe b zu Beginn 

25 des Iterationsprozesses zunachst mit Hilfe eines einfachen 

Matrix- Vektor-Produkts bestimmt werden. Die Matrix A wird 

aus Griinden der Effizienz vorzugsweise nicht explizit 

berechnet . Da samtliche Matrizen bei den CG-Verf ahren 

lediglich in Form von Matrix-Vektor-Produkten auftreten, ist 

30 es giinstiger, derartige Operationen sequentiell, d.h. in Form 

von zwei separaten Matrix-Vektor-Produkten, auszuwerten. 

-1 

Daruber hmaus geniigt anstelle der Inversen M die 
Berechnung einer geeigneten Faktorisierung von M. Die 
betroffenen Matrix-Vektor-Produkte miissen dann allerdings 
35 durch passende Riicksubs titut ions -Rout inen ersetzt werden. 
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Analog zur Vorkonditionierung von links wird die 
Vorkonditionierung von rechts durchgef uhrt : 



5 



A M 



-1. 



M x 



= b <=> A x 



= b 



mit A = 



A M 



-1 



x = M x 



(59) . 



Bei Anwendung der Gleichung (59) mufi die gesuchte Losung x 
nach dem eigentlichen Integrationsprozefc mit Hilfe eines 

10 zusatzlichen Matrix- Vektor-Produkts gemaS x = M _1 x berech.net 
werden . 



Da das Konvergenzverhalten des Iterationsprozesses durch die 
Verteilung der Eigenwerte von A bestimmt wird, sind beide 

15 Varianten der Vorkonditionierung bezuglich des 

Konvergenzverhaltens Equivalent. Im Vergleich zur links- 
Transformation hat die Vorkonditionierung von rechts 
allerdings den Vorteil, daS das Residuum aus Gleichung (54) 
als FehlermaS unverandert ubernommen werden kann. Im Fall der 

20 Vorkonditionierung von links la£t sich das Residuum aus 
Gleichung (54) meist nicht ohne Zusat zauf wand bestimmen. 
Vorteilhaft wird dann anstelle des in Gleichung (54) 
festgelegten Residuums die Form 



verwendet . Hier entstehen durch die Matrix M Verzerrungen, 
die das Abbruchkriterium ungiinstig beeinflussen konnen. 

30 JACOBI -VORKONDITIONIERUNG 




25 



rH = b — A xH = M^rH 



(60) 




-1 



Eine mogliche. Vorkonditionierung besteht aus der 
Diagonalmatrix 
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fAii fiir i = j 

M 



'ID 



Aaa fur i = j 

J (61) . 



sonst 



Die Anwendung dieser sogenannten Jacobi-Vorkonditionierung 
von links bzw. von rechts entspricht einer diagonalen 
Skalierung von A, bei der die Diagonal element e durch zeilen- 
bzw. spaltenweises Durchmultiplizieren des Gleichungssystems 
auf eins gebracht werden. 



10 Blockvariante der Jacobi-Vorkonditionierung. Hierbei muS eine 



Eine weitere vorteilhafte Vorkonditionierung ist eine 
Blockvariante der Jacobi-Vorkc 
Indexmenge der Freiheitsgrade 

F = {1,2, . . . ,N} 

15 zunachst gemafi 

F = U F v 



•w 

w 



in paarweise disjunkte Teilmengen F w aufgeteilt werden. Die 

20 Anzahl dieser Teilmengen sei mit W bezeichnet. AnschlieSend 
ergibt sich die Matrix M iiber die Beziehung 



[Ai j fiir i, j e F w 

Mij = \ n J (62) 
J 0 sonst 



25 Bei einer derartigen Block- Jacobi-Vorkonditionierung sind die 

Matrizen M und M Blockdiagonalmatrizen . Die Berechnung des 

-1 

Matrix-Vektor-Produkts M x laSt sich mit Hilfe separater 

Faktorisierungen der einzelnen Diagonalblocke effizient 
realisieren. 



Ein besonders gutes Konvergenzverhalten stellt sich ein, wenn 
die Indexmengen F w physikalischen Teilgebieten zugeordnet 
werden konnen. Die Erzeugung derartiger Teilgebiete, im 
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einfachsten Fall Wvirf el , ist durch schrittweises Unterteilen 
des zu untersuchenden Korpers in Form einer baumartigen 
Struktur moglich. Bei jedem Verf einerungsschritt werden die 
bereits generierten Wurfel vorzugsweise in bis zu acht 
5 Tochterwiirf el mit halber Seitenlange auf gespalten . Leere 
Gebiete werden ublicherweise ignoriert . Die grobste 
Verf einerungsstuf e, also ein einziger Block, wird als Stufe 0 
bezeichnet . Entsprechend erhoht sich der Stuf enindex bei 
jedem Verf einerungsschritt urn eins . 

10 

Vorteilhaft wird vor dem eigentlichen IterationsprozeS eine 
Permutation der Freiheitsgrade durchgefuhrt , damit die sich 
im Gleichungssystem ergebende Blockstruktur programmtechnisch 
handhabbar wird. Dazu werden Indizes derart umgeordnet, daS 
15 auf einanderf olgenden Freiheitsgraden gemeinsame Teilgebiete 
zugeordnet sind. Daraufhin konnen die den Indexmengen F w 
zugeordneten Teilproblem aus A extrahiert und mit Hilfe einer 

geeigneten Faktorisierung gelost werden. Die sich ergebenden 

. . -1 
Faktorisierungen werden dann zur Blockdiagonalmatrix M 

20 zusammengesetzt . 

In Fig. 6 wird beispielhaft die Einteilung einer 
Autokarosserie fur die Verf einerungsstuf en 2,3 und 4 gezeigt. 

25 Die Anwendung der Jacobi-Vorkonditionierung ist nicht auf 
physikalische Teilgebiete derselben GroSe beschrankt . 
Speziell bei einer Diskretisierung mit stark variierenden 
ElementgroSen bietet es sich an, den Korper in 
unterschiedlich groSe Teilgebiete zu unterteilen. In Fig. 7 

30 ist eine adaptive Block- Jacobi-Vorkonditionierung fur den 
Fall einer Leiterbahn uber einer unendlich ausgedehnten 
Massef lache dargestellt . Damit die sich ergebenden 
Stromverteilung durch die angesetzten Basisf unktionen 
moglichst genau nachgebildet werden kann, wird der Bereich um 

35 die Leiterschleif e fein vernetzt. Vorzugsweise wurden die 
Teilgebiete so gewahlt, daS in jedem Gebiet nicht mehr als 
300 Freiheitsgrade enthalten sind. 
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DAS FFS-VERFAHREN (FAST- FREQUENCY- STEPPING) 

Viele praktische EMV-Pr obi erne erfordern eine breitbandige 
Charakterisierung der zu untersuchenden Objekte. Sowohl 
5 Storf estigkeits- als auch Emissionsrichtlinien erstrecken 

sich iiber breite Frequenzbereiche , so daS zur Simulation der 
entsprechenden Messungen umf angreiche Rechenserien 
durchgefiihrt werden miissen, bei denen die Frequenz f, analog 
£ zur realen Messung, schrittweise von einer kleinsten zu 
"j^O untersuchenden Frequenz f m i n auf eine groSte zu untersuchende 
Frequenz f max erhoht wird. Damit kein kritischer 
Frequenzbereich iibersehen wird, sollte die verwendete 
Schrittweite Af moglichst klein gewahlt werden. 

15 Ein Ansatz zur Reduzierung des Rechenauf wands bei Verwendung 
iterativer Losungsverf ahren folgt aus der Beobachtung, daS 
sich die gesuchten Losungsvektoren fur genugend kleine 
Frequenzschritte nicht allzu sehr voneinander unterscheiden . 
Also kann der letzte der berechneten Losungsvektoren als 

2 0 neuer Startwert vorgegeben werden. 



Eine weitere vorteilhafte Auspragung besteht darin, den 
Startwert durch Extrapolation aus bereits berechneten 
Losungsvektoren zu bestimmen. Fur den Sonderfall 
25 aquidistanter Frequenzen laSt sich beispielsweise die 
Beziehung 

(0) 

x i = x i-3 _ 3x i-2 + 3x i-l (63) 



30 



zur quadratischen Extrapolation eines Startwerts aus den 

let z ten drei Losungsvektoren Xi_3/ x^_ 2 und x^.^ angeben. 



GR 97 P 1722 



34 

Ein schnelles Konvergenzverhalten kann beobachtet werden, 
wenn eines der Gleichungssysteme mit einem direkten 
Losungsverf ahren gelost und die berechnete Faktorisierung 
anschlieSend als Vorkonditionierer zur iterativen Losung des 
5 restlichen Gleichungssystems ubernommen wird. 

Durch Kombination dieses Schemas mit dem oben beschriebenen 
Extrapolationsansatz ergibt sich ein wirkungsvolles Vorgehen 
zur breitbandigen Analyse von elektrodynamischen 
10 Streuproblemen . Dieses Verf ahren wird im folgenden als FFS- 
Verf ahren bezeichnet . Der Algorithmus fur den Sonderf all 
aquidistanter Frequenzen lautet : 

f : = f min 

15 berechne A und b 

-1 

berechne M : = A (bzw. eine Faktorisierung von A) 
-1 

X : = M b 

solange f < f max 

f := f + Af 

2 0 berechne A und b 

extrapoliere einen Startwert i(°) 

lose AM x = b iterativ 
-1 . 

X := M x 



25 



Ende (Alg-l) 



Die Rechenserie wird beispielsweise an der tiefsten Frequenz 
f min gestartet. Nachdem das entsprechende Gleichungs system 
aufgestellt und die Matrix A geeignet faktorisiert wurde, 

30 ergibt sich der erste Losungsvektor durch Rucksubstitution . 
Die berechnete Faktorisierung von A wird anschlieSend als 
Vorkonditionierer M zur iterativen Losung der restlichen 
Gleichungssysteme ubernommen. Die Frequenz f wird dazu 
schrittweise von f m i n bis auf f max erhoht, wobei im Fall 

35 einer Vorkonditionierung von rechts beachtet werden muS, daS 
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nicht der eigentliche Losungsvektor x, sondern die 

transf ormierte GroSe x extrapoliert wird. Der eigentliche 

Losungsvektor x kann nach dem Iterationsprozefi uber die 

Beziehung x:=M -1 x bestimmt werden. 

5 

Der benotigte Rechenauf wand wird beim FFS-Verf ahren im . 

wesentlichen durch die Faktorisierung der Matrix M, d.h. 

durch das direkte Losen des Gleichungssystems , bestimmt. Mit 

-1 

Hilfe der Matrix M als Vorkonditionierer beschrankt sich 

10 die anschliefiende iterative Losung der restlichen 

Gleichungssysteme zumeist auf einige wenige Iterationen. 



Die Reihenfolge, mit der die einzelnen Frequenzen beim FFS- 
Verf ahren durchlaufen werden, mufi nicht notwendigerweise, wie 
15 in Algorithmus (Alg-1) dargestellt, von f m i n nach f max 

erfolgen. Ebenso sind Varianten des Verfahrens denkbar, die " 
an der hochsten Frequenz f max starten, das entsprechende 
Gleichungssystem direkt losen und anschliefiend die Frequenz 
schrittweise verringern. 



20 



it 



-l 

Weiterhin kann der Vorkonditionierer M fur eine mittlere 
Frequenz f c berechnet werden. Anschliefiend werden die 
Frequenzen wahlweise von f m i n nach f max , von f max nach f m i n 
Oder zunachst von f c nach f m i n und anschliefiend von f c nach 



25 fmax durchlaufen. Der Vorteil dieser Alternative besteht 



darin, dafi die breitbandige Beschleunigungswirkung der Matrix 
M in beide Frequenz richtungen ausgenutzt werden kann. 



MULTIPOLENTWICKLUNG DER RETARD I ERTEN POTENT I ALE 

30 Im Mittelpunkt einer schnellen Multipolmethode stehen 

Verfahren, mit denen sich die Potentiale und Felder einer 
gegebenen Stromverteilung naherungsweise durch einige wenige 
skalare Koef f izienten charakterisieren lassen. Wie durch den 
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Begriff "Multipolmethode" angedeutet, werden vorzugsweise 
spharische Multipolentwicklungen der Potentiale bzw. Felder 
eingesetzt . 

Da das in dieser Erfindung verwendete Linien- 

Gewichtungsverf ahren (ausf iihrliche Erlauterungen, siehe oben) 
auf den retardierten Potentialen A und cp aufsetzt, werden 
diese beiden GroSen durch entsprechende Multipolentwicklungen 
charakterisiert . 



Alternativ dazu ist es moglich, eine Multipolentwicklung der 
elektrischen Feldstarke E durchzufuhren. Im Vergleich zu den 
Multipolentwicklungen der retardierten Potentiale A und (p 
ergeben sich hierbei jedoch komplexere analytische 
15 Zusammenhange, die numerisch schwer zu handhaben sind. 

Gesucht ist eine Menge von Losungen der homogenen Helmholtz- 
Glei chung 



20 (v 2 + k 2 J • u = 0 



(64) 



in Kugelkoordinaten, wobei 
k die Wellenzahl, 

u eine Funktion, die die Helmhol tz -Gleichung erfullt, 

25 bezeichnen. 

Mit Hilfe dieser Losungen last sich jede Funktion u in Form 
einer unendlichen Reihe darstellen. Die dabei auftretenden 
Koef f izienten konnen haufig analytisch bestintmt werden. 

30 

Wird die Helmhol tz -Gleichung (64) mit Hilfe der in Fig. 8 
dargestellten Kugelkoordinaten r, S, a ausgedruckt, so folgt 
die Formulierung 

35 — — \rr — \+— sinS — + — + kTu = 0 

^drV drJ ^sxnS^V ddJ J sin 2 s ^2 
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(65) 



Zur Losung dieser partiellen Dif f erentialgleichung wird ein 
Separationsansatz 

5 

u = U r (r) • U$(S) - U a (oc) (66) 

gewahlt, bei dem die Faktoren U r , U$ und Ua jeweils nur von 
einer Koordinate abhangen. Nach Einsetzen von Gleichung (66) 
10 in Gleichung (65) und Multiplikation mit 

^ r 2 sin 2 5 

u 



15 



ergibt sich die Gleichung 

sin 2 S d ( 2 dU^ sin 3 d f . Q dU^ 1 d 2 ^ 



Uj- drV dry Us dS v dS -/ U a ^2 

(67) . 

Der a-abhangige Term laEt sich abspalten und mit Hilfe der 
2 0 Separationskonstanten m durch 

ft- 1 d 2 U a 2 , 

» a = (68) 



U a da 2 



ersetzen. Damit folgt aus Gleichung (67) die von a 
25 unabhangige Beziehung 

1 d ( 2 dU r ^ 1 d ( . _ dUa^ m 2 ,22 

r — + — sin 9 — -f- — + Y^sT = 0 

U r dr v dr/ sin8 d5 V dS / sin 2 s 

(69) . 

30 Zur Separation der S-Abhangigkeit werden die entsprechenden 
Terme in Gleichung (69) durch die Separationskonstante 1 
ausgedruckt : 
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U$ sin 3 dS 



^ dS J 



m 



sin 2 S 



« -1(1 + 1) 



(70) 



Ubrig bleibt die von r abhangige Dif f erentialgleichung 



1 d ( 2 dUr-^ , , 

— ^ - l(l + l) + k 

U r dr V dr ) v ' 



2 2 
r = 0 



(71) 



Mit Hilfe des Separationsansatzes aus Gleichung (66) wird die 
Helmholtz-Gleichung demnach in drei gewohnliche 
10 Dif f erentialgleichungen 



^( r2 ^) +[(kr)2 . 1(1 + 1) ]. Dr . 0 



(72a) 



15 



_J_ ( sin 3 ^ 
sin 3 d8 V d9 7 



1(1 + 1) - 



m 



sin 2 9 



XJ$ = 0 



d^U, 



da 



a. 2 

f- + m^U a = 0 



(72b) 



(72c) 



iiberf iihrt . 

20 Die Gleichung (72c) ist eine gewohnliche 

Dif f erentialgleichung mit konstanten Koef f izienten und den 
Losungen 



25 



U a ,m ~ e 



ima 



(73) 



Aufgrund der Nebenbedingung U a (a+27t) =U a (a) kann m nur 
ganzzahlige Werte annehmen, d.h. meZ (Z als Menge der ganzen 
Zahlen) . 
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Die Gleichung (72b) ist eine spezielle Form der Legendre- 
Dif f erentialgleichung . Ihre Losungen sind die zugeordneten 
Legendre-Funktionen erster und zweiter Art 



US.lm ~ PT(cosd), Qf(cosS) (74) . 

Nahezu alle Legendre-Funktionen weisen fur S=0 und S=7T 
Singularitaten auf und werden daher vorzugsweise 

vernachlassigt . Eine Ausnahme bilden die Funktionen P^fcos S) 
mit l ( meN 0 , die mit Legendre-Polynomen 1-ter Ordnung P^ iiber 
die Beziehung 



^-M-fi-^r^-^piw (75, 

v ' dx m 

verkniipft sind. Damit folgt 

pf(cos S) = ^(cos S) (76) 



sowxe 



Pffcos S) = 0 fur m>l (77) . 

Die zugeordneten Legendre-Funktionen erster Art von nullter 
bis zweiter Ordnung lauten: 

P,?(cos 8) = 1 B?(cos S) = - cos 2 S - - 

z 2 2 

P°(cos (*) = cos 3 e£(cos S) = -3 sin S cos £ 

Pj^cos S) = - sin 3 pf (cos S) = 3 sin 2 S 



Die Gleichung (72a) ist ein Besselsche Diff erentialgleichung 
mit den spharischen Bessel- Funktionen 
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U r ,l ~ jl(kr), yi(kr) (78) 

als Losung. Diese sind mit den bekannten zylindrischen 
5 Bessel-Funktionen J±{x) , (x) uber die Beziehungen 

jlto = ^ • Jl + 0,5W , Yl(x) = ^ • Y 1+0#5 (x) (79) 

verkniipft und lassen sich durch elementare Funktionen 
10 darstellen. So gilt beispielsweise 

. . . sin(kr) , , cos(kr) 

Do(kr) = n Yo( kr ) = , — 

kr kr 

. , x sin(kr) cos(kr) . , cos(kr) sin(kr) 

Di( kr ) = \ -^- L yi kr) = - \ ' ^—J- . 

(kr) 2 ( kr )2 kr 



15 



20 



Durch Kombination der spharischen Bessel-Funktionen ergeben 
sich die spharischen Hankel- Funktionen erster und zweiter Art 

h^(kr) = ji (kr) + i • yi(kr) (80a) 
h\ '(kr) = jx(kr) - i • yi(kr) (80b) , 



die beziiglich der Darstellung der Gesamtlosung u zu den 
Funktionen j^(kr), (kr) Equivalent sind. Je nach Art der 

2 5 Funktion u konnen vorteilhaft spharische Bessel-Funktionen 

bzw. spharische Hankel -Funktionen vorzuziehen. Bleibt die 
Funktion u am Ursprung beispielsweise beschrankt, so geniigt 
es, lediglich die Funktionen j^kr) anzusetzen, da die 
Funktionen yj (kr) und damit die spharischen Hankel -Funktionen 

3 0 hjj^kr) , h( 2 ^(kr) fur r->0 singular werden. 
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Es ergeben sich die gesuchten Losungen n^ m der homogenen 
Helmholtz-Gleichung in Kugelkoordinaten zu 

uim = b!(kr) • pH(cos S) - e ima (81) 

5 

wobei (kr) eine Linearkotnbination aus spharischen Bessel- 
bzw. Hankel-Funktionen 1-ter Ordnung bezeichnet und leN Q und 
meZ . 

10 Fur die Gesamtlosung u ergibt sich unter Berucksichtigung von 
Gleichung (77) die Reihenentwicklung 

oo l 

u = Z Z c ml • u m i = 
1 = 0 m = -l 

(82) 

= Z Z Cml • bi(kr) • pW(cos S) • e ima 
1 = 0 m = -l 

15 mit den konstanten Koef f izienten Ci m . 

^ KUGE L FLACHENFUNKT I ONEN 

Werden die winkelabhangigen Terme pj^(cos S) und e^" ma aus 

Gleichung (82) zusammengef afit , so ergeben sich die 

20 sogenannten Kugelf lachenf unktionen Y-J^S, a) . Mit einem 

geeigneten Normierungsf aktor versehen, folgt deren Definition 
zu 



25 



Y™(a,a) : = f" m j' . if (cob B) ■ e imx (83) 

x \ 4n (l + m) ! L 



fur 0 < m < 1. Kugelf lachenf unktionen mit m<0 lassen sich 
iiber eine Symmetriebeziehung 
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Yf m (S, a) = (-l) m Yf (S, a)* (84) 
definieren. Damit gilt eine Orthonomalitatsbedingung 

§ Y™ (S, a) ■ Y™'(S, a)* da = S ir • 8^ (85) 

a 

wobei 

H eine Oberflache einer Einheitskugei, 

8 ein Krone eke r- Symbol mit der Definition 

fl fur 1 = 1' 

5 11' = L (86) 
[0 sonst 

bezeichnen . 

Die Kugelf lachenf unktionen bilden auf CI einen vollstandigen 
Satz orthonormaler Funktionen bzgl . der Indizes 1 und m. Eine 
der wichtigsten Eigenschaf ten der Kugelf lachenf unktionen ist 
die daraus folgende Tatsache, daS sich jede beschrankte 
Funktion g(9,a) in eine Reihe nach den Kugelf lachenf unktionen 

Y m (S, a) entwickeln laSt: 

9(®.«) - £ Z C ml • Y™(» f a) (87) . 

1 = 0 m = -l 

Die Koef f izienten Cj m werden, unter Beriicksichtigung der 
Gleichung (85) , bestimmt iiber die Beziehung 

Cim = § g(S. «) • Yf (&, a)* da (88) 

n 

Die Kugelf lachenfunktionen nullter bis zweiter Ordnung 
lauten : 
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r-2 



15 2 q -i2a 

sin S • e 



3 27t 



,-1 _ !5 



sin $ • cos $ • e 



— xa 



,-1 



— • sin S • e- 1: 

871 



Y° = 



fl.fi 

V47I V2 



cos 2 9 - -I 



Y° = J— • COS S 
- 1 - V47C 



yJ = -J— • sin S • cos S • e ia 



Y, 1 = -J — • sin S • e 
^ V 87C 



la 



,,2 I !5 . 2 o i2a 
Yr = J • sin 9 • e 

2 « 3 27t 



ADDITIONSTHEOREM PER GREENSCHEN FUNCTION 
Die Greens che Funktion 

e -ikR 

G(r, r') = mit R = llr - r'll (89) 

R 

15 

beschreibt bis auf einen konstanten Faktor das retardierte 
\ Potential einer schwingenden Punktladung, die sich an einem 

Ort r' befindet. Die Entwicklung der Greenschen Funktion in 
spharische Losungen der Helmholtz-Gleichung hat folgende 

2 0 Form: 



oo 1 / -v 

-47cik £ ZiiiO^') Y i( s '> «') ^ M Y i( s > °0* Vr > r ' 

1 = 0 m=-l 

-47cik £ £ h^W) Y m (S' ( a') j!(kr) Y™(S, a)* Vr < r' 
1 = 0 m=-l 

(90) . 
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MULT I POLENTW I C KLUNGEN 

Mit Hilfe der Gleichung (90) la£t sich die 
Multipolentwicklung der retardierten Potentiale A und cp 
herleiten. Dazu wird Gleichung (90) in die Gleichungen (14a) 
5 und (14b) eingesetzt und die Reihenfolge von Integration und 
Summation vertauscht . Eine derartige Umformung ist zulassig, 
da der Sonderfall r=r' bereits in Gleichung (90) 
ausgeschlossen wurde . Die Entwicklung des Vektorpotentials A 
erfolgt durch separates Anwenden von Gleichung (90) auf die 
10 drei kartesischen Komponenten A x , Ay und A z . 

Je nach Lage der Beobachter- und Quellpunkte ergeben sich 
unterschiedliche Arten von Multipolentwicklungen. Die in 
Fig. 9a skizzierte Anordnung mit r>r' fuhrt zur sogenannten 
15 globalen Multipolentwicklung, bei der die Wirkung einer 
lokalen Quellverteilung G nach auSen in Form einer 
Multipolentwicklung dargestellt wird. Die beiden skizzierten 
Kreise mit den Radien d und D geben die Grenzen fur ein 
Nahfeld (kleiner als Radius d) und ein Fernfeld (groSer als 

2 0 Radius D) an. 

Die entsprechenden Reihenentwicklungen fur die retardierten 
Potentiale ergeben sich zu 

25 A(r) = -ilqio £ £ a? m h^kr) *f (», a)* 

1 = 0 m = -l 

(91a) , 

<p(-) = - 7 s £ £ p? m 4 2) m *±(*> «)* 

80 1 = 0 m = -l 

wobei sich die globalen Multipolkoef f izienten und p^ m 

3 0 uber die Beziehungen 

a lm = III J ( r ') 3l(kr') .Y"(S\ a') dV 
G 
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(91b) 

P?m = JJJ P(r') Dl(kr') YfCS', a') dV 
G 

berechnen lassen. Zur einfacheren Schreibweise wurden die 
drei skalaren Multipolkoef f izienten von A zu vektoriellen 
Koef f izienten zusammengef aSt . 

Der in Fig. 9b skizzierte Fall, bei dem sich die 
Beobachterpunkte naher am Koordinatenur sprung befinden als 
die Quellpunkte, d.h. r<r', fiihrt zu einer lokalen 
Multipolentwicklung . Hierbei wird die Wirkung einer gemigend 
weit entf ernten Quellverteilung G lokal entwickelt . Fur die 
retardierten Potentiale A und <p folgt 

A(r) = -i^ 0 Z £ j!(kr) Y»(», a)* 
1 = 0 m = -l 

(92a) 

cp(r) = - ±* £ 2 Pim Y i m ( S ' «)* 

80 1 = 0 m = -l. 

mit den lokalen Multipolkoef f izienten 

a L = JJJ J ( r ') 4 2) ( kr ') «') dV 

G 

(92b) . 

Plm = JJJ P( r ') 4 2) ( kr ') Y i m ( S '' a ') dV ' 
G 

Im Rahmen einer programmtechnischen Realisierung werden die 

Reihenentwicklungen (91a) und (92a) nach einer endlichen 

Anzahl von Termen abgebrochen. Gilt fur den auEeren 

Summations index 1=0 , 1 , . . . , L, so kann durch vollstandige 

Induktion gezeigt werden, daS die sich ergebenden 

2 

Multipolentwicklungen L-ter Ordnung stets (L+l) Terme 
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umfassen. Zur Multipolentwicklung der retardierten Potentiale 

2 

A, cp bis zur L-ten Ordnung werden somit insgesamt 4 (L+l) 
skalare Koef f izienten benotigt. 

Entsprechend dem Wert von 1 werden die wichtigsten 
Multipolkoef f izienten ublicherweise als 

1=0 : Monopolkoeff izient 

1=1: Dipolkoeff izient 

1=2 : Quadrupolkoef f izient 

1=3 : Oktupolkoef f izient 

1=4: Hexadekupolkoeff izient , usw. 

bezeichnet . 

Gemafi Fig. 9a und Fig. 9b wird das Beobachtergebiet von dem 
Quellgebiet durch zwei Kugelf lachen mit den Radien d bzw. D 
abgegrenzt. Fur d<<A. orientiert sich die Konvergenz der 
Multipolentwicklungen am statischen Fall und hangt im 
wesentlichen vom Verhaltnis d/D ab. Die Reihen konvergieren 
urn so besser, je kleiner d/D gewahlt wird. Im Fall elektrisch 
groSer Quellgebiete , d.h. d liegt im Bereich einer 
Wellenlange oder dariiber, bilden sich hingegen komplexe 
Interf erenzmuster zu deren Erfassung vorzugsweise L > kd 
gewahlt wird. 



MULTIPOLENTWICKLUNGEN DER TEILGEBIETE 

Im Rahmen einer schnellen Multipolmethode wird auf die oben 
beschriebene Unterteilung der zu untersuchenden Korper 
zuruckgegrif f en . 

In jedem aus der Unterteilung hervorgehenden Teilgebiet wird 
eine eigenstandige globale und lokale Multipolentwicklung 
angesetzt, d.h. fur jede Indexmenge der Freiheitsgrade F w 
exist iert sowohl eine Reihe globaler Multipolkoef f izienten 
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a lm ' Pirn als auch sine Reihe entsprechender lokaler 

Multipolkoef f izienten a^ m , p^ m . Die Multipolentwicklungen 

werden dabei beziiglich des jeweiligen Wiirf elmittelpunkts r w 
durchgef iihrt . 

5 

wird die Darstellung der Stromverteilung mit Hilfe der 
Basisfunktion vy n (siehe die Gleichungen (29a) und (29b) ) in 
die Gleichungen (91a, b) und (92a, b) eingesetzt, so folgt fur 
die globalen Multipolentwicklungen in den einzelnen 
10 Teilgebieten 



15 



A w (r) = -ikup J 2 a.\ m w 4 2 )(kr) Y»(«. a)* 



1 = 0 m = -l 



(93a) 



*»w = ~ dr 2 2 p ?m,w 4V) Y i m ( 9 ' «) 



^0 1 = 0 m = -l 



mit den globalen Multipolkoef f izienten 

a ?m, w = E ^ JJ Vn(r') Dl(kf') yf (S' f «') da' 
neF w Sn 



(93b) 



^ Plm, w = S Jnff V S ' • v n (r') ji(k?') Y™(§', a') da' 

neF w s n 

und fur die lokalen Multipolentwicklungen in den einzelnen 
Teilgebieten 

LI + 

25 A w (r) = -ikn 0 E 2>lm,w «) 

1 = 0 m = -l 

(94a) 

m ° 1 = 0 ra = -l 
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mit den lokalen Multipolkoef f izienten 




lm, w 



neF w \N W 



•n 



JJ n/ n (r') 4 2) {ki") Y™ (§', a') da' 



(94b) . 



Pirn, w 



n 



JJ V S ' • M> n (r') h^(kf') Y™(5', a') da' 



neF w \N W 



Die Tilde iiber den Kugelkoordinaten r, S und a deutet an, daS 
es sich hierbei urn lokale Koordinaten beziiglich des 

jeweiligen Wiirf elmittelpunkts handelt . Bei Bestimmung der 
lokalen Multipolkoef f izienten in Gleichung (94b) wurden 
Quellgebiete beriicksichtigt , die aufierhalb der Nachbargebiete 
liegen, d.h. neF\N w . 

Da die zu untersuchenden Korper im Bereich der EMV- Simulation 
(Leiterplatten, Verkabelungen, etc.) vorzugsweise auf kleine 
Teilgebiete fiihren, wird das Konvergenzverhalten der 
Reihenentwicklungen aus den Gleichungen (93a) und (94a) 
vorwiegend durch das Verhaltnis d/D bestimmt . Fur den Fall, 
daE die globalen Multipolentwicklungen (93a) auSerhalb der 
Nachbargebiete ausgewertet werden, gilt bei wiirf elf ormigen 
Teilgebieten in drei Raumdimensionen sowohl fur die globalen 
als auch fur die lokalen Entwicklungen 

d 1 
D Vi " 



STAB I L I TATS PROBLEME BEI TIEFEN FREOUENZEN 

Erfolgt die Einteilung des zu untersuchenden Korpers nach 
Freiheitsgraden bzw. nach Basisf unktionen, so konnen bei 
tiefen Frequenzen Instabilitaten beziiglich der globalen 
Multipolentwicklung des elektrischen Skalarpotentials <p 
auf treten . 
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Fia. 10 verdeutlicht die Problematik am Beispiel einer 
einfachen Leiterschleif e , deren Skalarpotential durch vier 
globale Multipolentwicklungen dargestellt werden soil 
5 (Fig.lOa) . Mit fallender Frequenz verschwindet die 
Gesamtlandung auf der Leiteroberf lache, da sich eine 
zunehmend gleichmaEigere Stromverteilung einstellen wird. 
Gleichzeitig wachsen jedoch aufgrund der 
Kontinuitatsgleichung die Ladungen in den einzelnen 

10 Basis funktionen mit l/o (vergleiche Gleichung (29b) ) . Solange 
sich die Basisf unktionen iiberlappen, konnen sich die Ladungen 
gegenseitig ausldschen und die Forderung nach einer 
verschwindenden Gesamtladung erfiillen. Durch das 
"Auseinanderziehen" der Basisf unktionen bleiben jedoch auf 

15 den Grenzelementen Ladungen bestehen, die sich nicht mehr 
gegenseitig kompensieren konnen, da sie unterschiedlichen 
Multipolen zugeordnet sind (siehe Fig.lOb). Wird die 
Multipolentwicklung von cp nach endlich vielen Termen 
abgebrochen, so ergeben sich Abbruchf ehler , die mit fallender 

2 0 Frequenz zunehmend Probleme verursachen. 

Abhilfe wird dadurch geschaffen, daS den Ladungen Gelegenheit 
zur Ausloschung gegeben wird, bevor die eigentlichen 
Multipolentwicklungen durchgefiihrt werden . Eine mogliche 
jlQ5> programmtechnische Realisierung besteht darin, samtliche 
globale Multipolentwicklungen nach Elementen anstatt nach 
Basisf unktionen durchzuf uhren . 



DIE SCHNELLE MULTI POLMETHODE (FMM) 

30 Die schnelle Multipolmethode , nachfolgend als FMM (engl.: 

fast multiple method) bezeichnet, ist z.B. aus [4] bekannt . 
Samtliche Implement ierungen der FMM fur den dynamischen Fall 
eignen sich fur Streuprobleme, bei denen Teilgebiete 
wesentlich groSer sind als die Wellenlange. Bei derartigen 

35 Problemen ist eine hohe Ordnung der verwendeten 
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Multipolentwicklungen notwendig. Ein wesentlicher Nachteil 
bei diesetn Formalismus ergibt sich aus der Eigenschaf t, dafi 
die Multipolkoef f izienten nicht explizit berech.net werden. 
Anstelle der expliziten Berechnung wird die zu entwickelnde 
Funktion an diskreten Punkten auf einer umgebenden 
Kugelflache "abgetastet " , woraufhin die Testwerte unmittelbar 
weiterverarbeitet werden. Da eine Unterdriickung der 
Multipolterme hoherer Ordnung dabei nicht ohne wei teres 
moglich ist, konnen storende Aliasing-Ef f ekte auftreten, 
welche unter Umstanden zu Fehlern in den 
Multipolentwicklungen fuhren. 

Das hier vorgestellte Verfahren eignet sich insbesondere zur 
Simulation von Streuproblemen, bei denen die geometrischen 
Strukturen der Modelle und damit auch die Dimensionen der 
Teilgebiete kleiner als die Wellenlange sind. Unter 
derartigen Voraussetzungen bleiben die Interf erenzef f ekte 
begrenzt, so dafi bereits wenige Multipolkoef f izienten zur 
Beschreibung von Potent ialen bzw. Feldstarken ausreichen. Mit 
dem hier beschriebenen Stabilisierungsverf ahren verhalt sich 
das Verfahren auch fur tief e Frequenzen stabil . 

Besonders interessant ist der Einsatz der schnellen 
Mult ipolmethode im Bereich von strahlungsgebundenen EMV- 
Problemen. Die zu simulierenden Messungen werden vorzugsweise 
auf einige dB genau durchgef iihrt . 

Bei Anwendung der beschriebenen iterativen Losungsverf ahren 

auf voll besetzte Gleichungssysteme ergibt sich durch das zu 

berechnende Matrix- Vektor-Produkt vorzugsweise eine 

2 

Komplexitat von 0(N ) Fliefikommaoperationen pro 

Iterationsschritt . Der Einsatz der FMM ermoglicht hingegen, 

aufgrund der besonderen Struktur der zugrundeliegenden 

Integralgleichungen das Matrix- Vektor-Produkt und damit auch 

2 

einen Iterationsschritt in weniger als 0(N ) 
Fliefikommaoperationen durchzuf iihren . 
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Ausgegangen wird von einer komprimierten Darstellung der 
Matrix Z in der Form 

Z = Z' + LTG (95) , 

5 

wobei 

Z' denjenigen Teil der Matrix Z, der die Kopplung 

zwischen benachbarten Teilgebieten beschreibt, 
LTG den restlichen Teil der Matrix Z, der die Kopplung 

10 zwischen weit entfernten Teilgebieten beschreibt, 

d.h.ieF w ,jeF\N w 
^ bezeichnen. 

Fiir den Fall, dafi auf die oben erlauterten 
15 Stabilisierungsmafinahmen (siehe Kapitel : Stabilitatsprobleme 
bei tiefen Frequenzen) verzichtet wird, bestimmt sich die 
Matrix Z' gemaS 

fiir i e F w , j e N w 

(96) . 

sonst 



20 




Im stabilisierten Fall treten zusatzliche Korrekturterme auf 



In Fig. 11 ist ein Beobachtergebiet BG der Indexmenge F w 
dargestellt, das von Nachbargebieten NG der Indexmenge N w 
25 umgeben ist. Die Nachbargebiete NG sind somit direkt 

gekoppelt. Weit entfernte Gebiete WEG sind als indirekte 
Kopplungen zum Beobachtergebiet BG dargestellt . Globale 
Multipolentwicklungen GMPE werden in Gebieten mit einem Kreuz 
durchgef tihrt , wohingegen die lokale Multipolentwicklung in 
30 dem Beobachtergebiet BG durchgef iihrt wird. 

Nachfolgend werden die globalen Multipolkoef f izienten a^ m , 

a g 
Plm w der einzelnen Teilgebiete in einem Vektor c 

1 

zusammengef afit . Analog dazu umf aSt ein Vektor c alle lokalen 
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Multipolkoef f izienten a^ m w , p^ m w . Die Nachbargebiete in 

Fig. 11 bleiben im Rahmen der lokalen Multipolentwicklung 
unberucksichtigt , da die entsprechenden Beitrage bereits in 
der Matrix Z' enthalten sind. Werden die 
5 Multipolentwicklungen bis zur L-ten Ordnung durchgef uhrt , so 
gilt im Fall von W Teilgebieten 

c9 c 1 e c 4W ( L+1 ) 2 . 
10 Die Vektoren c g und c"'' berechnen sich iiber die Beziehungen 
c g = G • I (97a) 
c 1 = T • c g = T • G • I (97b) , 

15 

d.h. die Matrix G ermoglicht die Bestitnmung der globalen 
Multipolentwicklungen in den Teilgebieten fur eine gegebene 
Stromverteilung I. AnschlieEend berechnet die sogenannte 
Trans 1 at ionsmatrix T daraus die lokalen Multipolkoef f izienten 
20 c . Die Translationsmatrix T sammelt die in Fig. 11 durch 

Kreuze dargestellten globalen Multipolentwicklungen in einer 
lokalen Multipolentwicklung. SchlieSlich werden die lokalen 
Multipolentwicklungen mit Hilfe der Matrix L an den 
Beobachterpunkten ausgewertet und zu den Nachbarbeitragen 
25 addiert : 

ZI = Z'l + Lc 1 (98) . 

Im nicht stabilisierten Fall sind die Matrizen G und L 
30 Blockdiagonalmatrizen, wobei die einzelnen Diagonalblocke die 
Multipolentwicklungen der Teilgebiete mit den darin 
enthaltenen Basisfunktionen verkniipfen. Beziehungen zur 
Berechnung der Matrixelemente von G folgen unmittelbar aus 
der globalen Multipolentwicklung gemafi den Gleichungen 
35 (93a, b) . Die gesuchten Matrixelemente entsprechen den in 
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Gleichung (93b) auf tretenden Integral termen . Beziehungen zur 
Berechnung der Mat rixel entente von L ergeben sich durch 
Einsetzen der Gleichung (94a) in die Gewichtungsgleichung 
(34) . Hierbei wird vorzugsweise beachtet, daS die GroSen A, cp 
5 aus Gleichung (34) durch die Beitrage der jeweiligen 
Teilgebiete A^, (p w zu ersetzen sind. 

Impedanzbehaf tete Elemente lassen sich gemaB der oben 
beschriebenen Vorgehensweise beriicksichtigen . Da sich 

10 lediglich fur benachbarte Elemente Korrekturterme ergeben, 
bleiben durchzuf iihrende Modif ikationen auf die Matrix Z' 

j^j beschrankt . 

Die Berechnung der Translationsmatrix T wird nachfolgend 
15 beschrieben. 



TRANS LATIONSOPERATOREN 

Zur eindeutigen Charakterisierung einer Multipolentwicklung 
gehort neben den einzelnen Mulipolkoef f izienten die 
20 Spezif ikation des Koordinatensystems , beziiglich dessen die 
Entwicklung durchgefuhrt worden ist . Die 

Multipolkoeff izienten sind dabei unmittelbar sowohl vom 
W y Entwi cklungs z ent rum (Koordinatenursprung) als auch von der 

Definition der WinkelgroSen & und a (Orientierung der 
25 Koordinatenachsen) abhangig. 

Die im Rahmen der schnellen Multipolmethode benotigten 
Opera toren zur Umrechnung der Koeff izienten einer gegebenen 
Multipolentwicklung auf ein neues Koordinatensystem sind 
3 0 Translationsoperatoren, da die Orientierung der 

Koordinatenachsen, im Gegensatz zu den sogenannten 
Rotationsoperatoren, erhalten bleibt. Falls die 
Multipolentwicklungen nach endlich vielen Koeff izienten' 
abgebrochen werden, lassen sich die gesuchten 
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Translationsoperatoren in Form voll besetzter quadrat ischer 
Translationsmatrizen angeben. 



Die nachfolgend verwendeten geometrischen Grofien sind zur 
5 Veranschauli chung in Fig. 12 dargestellt. Der Ortsvektor eines 
Punktes P wird im alten Koordinatensystem mit f = (r, 5, a) und 
im neuen Koordinatensystem mit r = (r, S, a) bezeichnet . 
Weiterhin wird der alte Koordinatenursprung 0 durch den 
Ortsvektor r' = (r', S', a') im neuen Koordinatensystem eindeutig 
10 f estgelegt . Die Translation erf olgt damit in negativer r' - 
Richtung. Hier werden vorzugsweise solche 

Translationsoperatoren betrachtet, die eine gegebene globale jQff 
Multipolentwicklung in eine neue lokale Multipolentwicklung 
umrechnen, also r < r' . 

15 

Ausgegangen wird bei der Herleitung der gesuchten 
Translationsoperatoren von dem Additionstheorem der 
Greenschen Funktion (siehe Gleichung (90) ) , woraus sich die 
Beziehung 

20 

h^kr) Y°(3, a) = ^ J) 2 hfV) Yf(S', a') ^(kr) Y^S, a)* 

1=0 m=-l 

(99) 

zur Translation von globalen Monopoltermen (l=m=0) ergibt 

2 5 Hierbei sind die Identitaten 

hfV) - i ^ und Y°(S, a) . 

berucksichtigt worden. Die Tatsache, daS Gleichung (99) die 

3 0 gewunschte Translation beschreibt, laSt sich daran zeigen, 

dafi auf der linken Seite der alte Monopolterm steht, der 
durch eine Uberlagerung aus neuen lokalen Multipoltermen 
ersetzt werden kann. Die einzelnen Vorfaktoren 
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4^ h^(kr') Y™(S' ( a') 



hangen nur vom Translationsvektor r' ab und konnen im voraus 
berechnet und in der ersten Spalte der Translationsmatrix . 
5 abgelegt werden. 

Fur eine programmtechnische Realisierung werden vorzugsweise 
die gesuchten Operatoren rekursiv durch sukzessive 
Differentiation der Gleichung (99) bestimmt. 



10 



Zur Beschreibung der benotigten Dif f erentiationsbeziehungen 
werden zunachst die Dif f erentialoperatoren 



15 



d+ = 



ax 



+ i 



5x 



- i 



d 

dy 

_d_ 
dy 



(100a) 



(100b) 



a. = 



dz 



(100c) 



20 eingefiihrt. Damit gilt beispielsweise 



25 



hjftkr) Y 0 °(9, a)l = £ k 4 2 )(kr) Y^S, a) 



hjftkr) Y°(S, a) 
h^(kr) Y 0 °(S, a) 



= - 1 J|k4 2 )(kr)Y l - l (S,a) 
= k h^kr) Y°(S, a) 



(101a) 



(101b) 



(10lc) , 



30 



d.h. Multipolterme erster Ordnung (Dipole) lassen sich direkt 
durch geeignetes Dif f erenzieren des Monopol terms herleiten. 
Entsprechende Beziehungen fur Multipolterme hoherer Ordnung 
werden im nachfolgend hergeleitet. 
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Fur spharische Hankel-Funktionen lassen sich die elementaren 
Dif f erentiationsbeziehungen 



= kr 1+1 4 2 l 1 (kr) sin S e ict 



(102a) 



d- 



= kr 1+1 4 2 2 1 (kr)sinSe- ict 



(102b) 



r 1 + 1 4 2) (kr)l = kr 1+1 hfl^kr) cos S 



(102c) 



10 r.v 

zeigen. Entsprechendes gilt fur Kugelf lachenf unktionen 



Yf (» f a) 



.1+1 



YfJi 1 ^. «) 



= a™(21 + l) 



L l + 1 



J. +2 



(103a) 



15 



Yf(3, a) 



vf+M a) 



= -b™(21 + l) 



a +2 



(103b) 



Yf (3, a) 



.1 + 1 



-C™(21 + l) 



Yf +1 (S, a) 



.1+2 



(103C) 



mit den konstanten Koef f izienten 



20 



m _ j (l + m + l)(l + m + 2) 
(21 + l)(21 + 3) 



a l = 



(104a) 



btn _ 
T = 



|(l - m + l)(l - m + 2) 



(21 + l)(21 + 3) 



(104b) 



25 



(l + m + l)(l - m + l) 
(21 + l)(21 + 3) 



(104c) , 
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die sich aufgrund des Normierungsf aktors aus der Definition 
der Kugelf lachenf unktionen (siehe Gleichung (83)) ergeben. 
Durch die Kombination von Gleichungen (102a-c) rait 
Gleichungen (103a-c) folgen die Dif f erentiationsbeziehungen 
fur Multipolterme hdherer Ordnung 



(2), 



5 + [4 2 )(kr) yf («, a)] = a™(21 + l) Y^(S, a) + 



+ khfl^kr) sin S e ia Y^(S, a) (105a) 



ia „m/ 



a_[hP(kr) Y™(&, a) 



= -b m (21 + 1) 1 ^ ' Y™ a) + 
+ khfl^kr) sin 8 e _ia Y™^, a) (105b) 



15 



20 



25 



h^(kr) Yf (S, a) 



= -c 



^(21 + 1) Y m +i(9/ a) 



+ 



+ kh^ 2 2 1 (kr) cos S Y^(S, a) 



(105c) , 



die sich mit Hilfe der Rekursionsbeziehung fur spharische 
Hankel - Funkt ionen 



4 2 2 l0 cr) + hf) l(kr) - (21 + l) ^ 



(106) 



sowie den Rekur s i onsbe z iehungen fur Kugelf lachenf unktionen 



sin 9 e lct Y™(8, a) = a^" 1 • Y^(S, a) - a m • Yf+^S, a) 



m+1/ 



(107a) 



sin S e _ia Yf (S, a) = -b^ 1 • vf." 1 ^, a) - bj 1 ■ Yf^S, a) 



(107b) 



cos S Y™(S, a) = cf^ • Yf.^S, a) - cf ■ Y? +1 (9, a) (107c) 
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weiter umformen lassen. Dabei ergeben sich schlieSlich die 
symmetrischen Different iationsbeziehungen 



7 + 



4 2 >(kr) Y™($, a)l = a™ k hfil yf + + 1 1 (» f a) + 



+ a. 



(108a) 



4 2) (kr) Yf (S, a) 



= -b m k hfix Yf a) - 
- b^ 1 k 4 2 liM YfJi 1 ^. a) 



(108b) 



3 z [4 2 )(kr) Y^S, a)] = -c m k 4 2 i 1 Y™ +1 (»,a) 



+ c5 n _ 1 k4 2 2 1 MYf_ 1 (S,a) 



(108c) 



mit deren Hilfe sich die gesuchten Translationsoperatoren 
berechnen lassen. Hierzu wird durch sukzessives 
Dif f erenzieren von Gleichung (99) bezuglich der Koordinaten 
r', a' von Translationsoperatoren (1-1) -ter und 1-ter 
Ordnung auf Translationsoperatoren (1+1) -ter Ordnung 
geschlossen. Dabei sei beachtet, da£ d = -d' gilt. 



MEHRSTUF IGE ALGOR I THMEN 

Neben der beschriebenen einstufigen FMM existieren diverse 
mehrstufige Variant en, bei denen Teilgebiete 
unterschiedlicher GroSe verwendet werden. Zur 
2 5 Veranschauli chung sind in Fig. 13 relevante Teilgebiete fur 
den Fall einer zweistufigen FMM-Variante dargestellt. 

In Fig. 13 ist ein Beobachtergebiet BG der Indexmenge Fw 
dargestellt, das von Nachbargebieten NG der Indexmenge N w 
30 umgeben ist. Die Nachbargebiete NG sind somit direkt 

gekoppelt. Weit entfernte Gebiete WEG sind als indirekte 
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Kopplungen zum Beobachtergebiet BG dargestellt. Globale 
Multipolentwicklungen GMPE werden in Gebieten mit einem Kreuz 
durchgef uhrt , wohingegen die lokale Multipolentwicklung in 
dem Beobachtergebiet BG durchgef iihrt wird. 

Die Grofie der Teilgebiete wird bei mehrstufigen FMM- 
Algorithmen vorzugsweise proportional zum Abstand zwischen 
Beobachtergebiet und Quellgebiet bestimmt . Die 
Multipolentwicklungen auf den verschiedenen 
Verf einerungsebenen lassen sich mit Hilfe von speziellen 
Translationsoperatoren ineinander umrechnen, so daS 
vorteilhaft nicht fur jedes Teilgebiet der hierarchischen 
Gebietsstruktur eine vollstandige Multipolentwicklung 
durchgef iihrt wird. 

Im dynamischen Fall werden vorzugsweise fur grofiere 
Teilgebiete Multipolentwicklungen hoherer Ordnung 
durchgef iihrt . 

In Fig.l ist ein Blockdiagramm dargestellt, das Schritte 
eines Verfahrens zur Bestimmung des elektromagnetischen 
Feldes umfaSt . In einem Schritt la wird die globale 
Multipolentwicklung fur den in Teilbereiche unterteilten zu 
untersuchenden Korper, wie oben beschrieben, durchgef iihrt . In 
einem Schritt lb wird die lokale Multipolentwicklung fur den 
in Teilbereiche unterteilten Korper durchgef iihrt . SchlieSlich 
wird in einem Schritt lc durch eine Uberlagerung aus der 
globalen und der lokalen Multipolentwicklung das 
elektromagnetische Feld des Korpers bestimmt. 
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Im Rahmen dieses Dokuments wurden folgende Verof f entlichungen 
zitiert : 

[1] Sadasiva M. Raoi, Donald R. Wilton, Allen W. Glisson: 

Electromagnetic Scattering by Surfaces of Arbitrary 

Shape . IEEE Trans . Antennas Propagat . , Vol . 30 , No . 3 , 
pp. 409-418, May 1982. 

[2] Ning Yan Zhu and Friedrich M. Landstofer: Application of 
Curved Parametric Triagular and Quadrilateral Edge 
Elements in the Moment Method Solution of the EFIE. IEEE 
Microwave and Guided Wave Letters, Vol.3, No . 9 , pp. 3 19- 
321, September 1993. 

[3] Youcef Saas and Martin H. Schultz : GMRES : A Generalized 
Minimal Residual Algorithm for Solving Nonsymmetric 
Linear Systems. SIAM J. Sci. Stat. Comp., Vol.7, No. 3, 
pp. 856-869, July 1986. 

[4] V. Rokhlin: Rapid Solution of Integral Equations of 
Classical Potential Theory. Journal of Computational 
Physics, Vol.60, pp. 187-207, 1985. 
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Patentanspruche 

1. Rechnergestiitztes Simulationsverf ahren zur Bestimmung 
eines elektromagnetischen Feldes eines Korpers, der 
mehrere Teilgebiete umfaSt und mehrere Ladungen und 
Strome enthalt, 

a) bei dem in jedem der mehreren Teilgebiete eine globale 
Multipolentwicklung, die die Wirkung der Ladungen und 
Strome fur entfernte Punkte in Form einer 
Multipolentwicklung darstellt, und eine lokale 
Multipolentwicklung, die die Wirkung der Ladungen und 
Strome an Punkten innerhalb des jeden der mehreren 
Teilgebiete in Form einer Multipolentwicklung 
darstellt, durchgefiihrt werden, 

b) bei dem durch eine Uberlagerung aus der globalen 
Multipolentwicklung und der lokalen 
Multipolentwicklung der mehreren Teilgebiete das 
elektromagnet ische Feld des Korpers bestimmt wird. 

2. Verf ahren nach Anspruch 1, 

bei dem bis zu einem vorgebbar kleinen FehlermaS folgende 
Schritte zur Bestimmung des elektromagnetischen Feldes 
des Korpers iterativ durchgefiihrt werden, wobei als eine 
Anf angsbedingung 1=0 gilt: 

a) Berechnung der globalen Multipolentwicklungen mit 
globalen Multipolkoef f izienten gemaS 

C 9 = GI , 

wobei 

ein Vektor ist, der die globalen 

Multipolkoef f izienten der mehreren 

Teilgebiete umfaSt, 
I ein Vektor ist, der eine gegebene 

Stromverteilung angibt, 
G eine Matrix ist, die eine Bestimmung der 

globalen Multipolkoef f izienten in dem einen 
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Teilgebiet der mehreren Teilgebiete fur die 
gegebene Stromverteilung I angibt; 
b) Berechnung der lokalen Multipolentwicklung mit lokalen 
Multipolkoef f izienten gemaS 



c 



1 = T C 9, 



wobei 

c 1 ein Vektor ist, der die lokalen 

10 Multipolkoef f izienten der mehreren 

Teilgebiete umfafit, 
T eine Translationsmatrix ist, durch die die jjj^f 

globalen Multipole in lokale Multipole 
zusammengef iihrt werden; 
15 c) Bestimmung des elektromagnetischen Felds aus 

ZI = Z'l + Lc 1 , 



wobei 

2 0 Z eine Impedanzmatrix, 

Z' einen Teil der Impedanzmatrix Z, der 

Kopplungen zwischen den Teilgebieten 
darstellt , 

L eine Matrix zur Auswertung der lokalen 

25 Multipolkoef f izienten, 

bezeichnen . 



m 



3. Verfahren nach einem der vorhergehenden Anspriiche, 
bei dem die Teilgebiete gleich grofi sind. 

30 

4 . Verfahren nach einem der Anspriiche 1 Oder 2 , 

wobei die Teilgebiete in ihrer GroSe proportional zu 
einer Entfernung von einem Beobachtergebiet sind. 

35 5. Verfahren nach einem der vorhergehenden Anspruche, 

bei dem das eine Teilgebiet der mehreren Teilgebiete 



jeweils einem Bereich mit einheitlicher physikalischer 
Eigenheit zugeordnet wird. 

Verfahren nach einem der vorhergehenden Anspruche, 

bei dem in einem Verf eine rungs schritt das Teilgebiet der 

mehreren Teilgebiete in bis zu acht Bereiche aufgespalten 

wird. 

Verfahren nach einem der Anspruche 2 bis 6, 
bei dem ein impedanzbehaf tetes Element, das ein 
Bestandteil des Teilgebiets des Korpers ist, direkt in 
der Matrix Z' als Impedanz beriicksichtigt wird. 

Verfahren nach einem der vorhergehenden Anspruche, 
bei dem das elektromagnetische Feld fur vorgebbare 
Frequenzen bestimmt wird. 

Verfahren nach Anspruch 8 , 

bei dem die vorgebbaren Frequenzen durch eine minimale 
Frequenz und durch eine maximale Frequenz bestimmt sind, 
wobei bei der minimalen Frequenz gestartet wird und bis 
zu der maximalen Frequenz fortgesetzt mit einer 
vorgebbaren Schrittweite jeweils das elektromagnetische 
Feld bestimmt wird. 

Verfahren nach Anspruch 8, 

bei dem die vorgebbaren Frequenzen durch eine minimale 
Frequenz und durch eine maximale Frequenz bestimmt sind, 
wobei bei der minimalen Frequenz gestartet wird und bis 
zu der maximalen Frequenz fortgesetzt mit einer 
vorgebbaren Schrittweite jeweils das elektromagnetische 
Feld bestimmt wird. 

Verfahren nach Anspruch 8, 

bei dem die vorgebbaren Frequenzen durch eine minimale 
Frequenz und durch eine maximale Frequenz bestimmt sind, 
wobei bei einer Frequenz zwischen der minimalen Frequenz 
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und der maximalen Frequenz gestartet wird und bis zu der 
maximalen Frequenz oder bis zu der minimalen Frequenz 
fortgesetzt mit einer vorgebbaren Schrittweite jeweils 
das elektromagnetische Feld bestimmt wird. 

5 

12. Verfahren nach einem der vorhergehenden Anspriiche, 
bei dem Stabilitat bei tiefen Frequenzen gewahrleistet 
wird, indem die globale Multipolentwicklungen 
elementweise anstatt nach Basisf unktionen durchgefiihrt 

10 werden . 

13. Verfahren nach einem der vorhergehenden Anspriiche, 
bei dem eine elektromagnetische Vertraglichkeit des 
Korpers bestimmt wird. 



15 
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Zusammenf as sung 

Rechnergestiitztes Simulationsverf ahren zur Bestimmung eines 
elektromagnetischen Feldes eines Korpers, der mehrere 
Teilgebiete umfaSt und mehrere Ladungen und Strome enthalt 

Das Verfahren ermoglicht die Bestimmung eines 

elektromagnetischen Feldes eines Korpers, indem dieser Korper 
in Teilbereiche unterteilt wird und fur die Teilbereiche 
Multipolentwicklungen durchgefuhrt werden, die iiberlagert das 
elektromagnetische (Gesamt- ) Feld ergeben. Die sich ergebenden 
Dif f erentialgleichungen fur die Teilbereiche werden in 
Matrixf orm angetragen und numerisch gelost . Anwendung ist 
u.a. die Bestimmung der elektromagnetischen Vertraglichkeit 
(EMV) eines Korpers . 



Figur 1 
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